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Представлена математическая модель, описывающая основные физико-химические процессы, происходящие в системе газовая 
фаза ~ водная фаза ~ грунт. Модель описывает многокомпонентную фильтрацию водной фазы, движение газовой фазы, испа- 
рение (конденсацию) и кристаллизацию (плавление), диффузию, сорбцию и десорбцию радионуклидов, теплопроводность и 
конвективный теплообмен. На основе модели было создано проблемно-ориентированное программное обеспечение, которое 
позволяет прогнозировать миграцию радионуклидов в поверхностном слое почвы. Приведены результаты моделирования вер- 
тикальной миграции радиоактивных веществ в поверхностном слое почвы. 


Изучение поведения радионуклидов в почвен- 
ном покрове является важной задачей с точки зре- 
ния экологической безопасности регионов, в кото- 
рых расположены предприятия атомной промы- 
шленности. Радионуклиды могут попасть в почву в 
результате внештатных ситуаций, связанных с ат- 
мосферными выбросами и проливами радиоактив- 
ных веществ. Кроме этого в 50-70 годы XX в. для 
хранения жидких радиоактивных отходов (ЖРО), 
образующихся в процессе работы предприятий 
ядерного топливно-энергетического цикла, ис- 
пользовались поверхностные хранилища предста- 
вляющие собой естественные (оз. Карачай) или 
специально созданные гидроизолированные водо- 
емы [1,2]. Миграция радионуклидов из естествен- 
ных и поверхностных водоемов представляет собой 
потенциальную угрозу для экосистемы. 

Миграция радионуклидов в почвенном покрове 
определяется значительным количеством взаимо- 
связанных физико-химических процессов, проис- 
ходящих в системе: газовая фаза - водный раствор 
- вмещающая порода. Вследствие сложности этих 
процессов для описания поведения радионуклидов 
целесообразно использовать методы математиче- 
ского моделирования. С помощью математическо- 
го моделирования возможно описать основные 
процессы, происходящие в системе, с учетом их 
нелинейного влияния друг на друга. 

В настоящей работе представлена математиче- 
ская модель, описывающая вертикальную мигра- 
цию радионуклидов в ненасыщенной зоне почвы с 
учетом фильтрации воды с поверхности почвы, 
фильтрации воздуха в почве, испарения (конденса- 
ции) и кристаллизации (плавления) воды, диффу- 
зии, сорбции и десорбции радионуклидов, радио- 
активного распада, энерговыделения в результате 
радиоактивного распада, конвективного и кондук- 
тивного теплообмена. Кроме этого, учитывается 
внешнее воздействие со стороны окружающей сре- 
ды (изменение атмосферного давления, температу- 
ры и влажности воздуха на поверхности почвы, а 
также атмосферные осадки). Приводятся результа- 
ты моделирования эволюции состояния поверх- 
ностного слоя почвы и вертикальной миграции ра- 
дионуклидов. 


В модели рассматривается одномерная верти- 
кальная миграция радионуклидов в поверхностном 
слое почвы, который представляет собой многофаз- 
ную систему жидкость - газ - твердая фаза. Рассма- 
триваемая система состоит из семи частей: вмещаю- 
щая порода (0), жидкая фаза (1), лед (2), газовая фа- 
за (3), граница раздела фаз (ГРФ) порода - жидкая 
фаза (4), ГРФ порода - лед (5) и ГРФ порода - газо- 
вая фаза (6). Для количественного определения до- 
ли порового пространства, занимаемого частью си- 
стемы Ф, используется понятие насыщенности 5 Ф : 

А Ф =-^-, (Ф = 1-3), 
т V 


где т - пористость; Ѵ ф - объем занимаемый фазой 
Ф; V - объем среды. Части системы с 1 по 3 зани- 
мают весь объем порового пространства: 

1=ІХ 

Ф=1 

Рассматривается поведение воды и радиону- 
клидов (индекс воды равен 1, радионуклиды нуме- 
руются с 2 по АО. Полагается, что радионуклиды не 
переходят в газовую фазу при испарении воды. При 
замерзании воды радионуклиды остаются на по- 
верхности вмещающей породы, контактирующей 
со льдом. Таким образом, радионуклиды могут на- 
ходиться только в жидкой фазе и в адсорбирован- 
ном состоянии на поверхности породы. В жидкой 
фазе количество /-ого радионуклида описывается 
объемной активностью й[: 

а[ = А = ^і_ 

1 Ѵ } тЗу 


где Ч'ф - активность /-ого радионуклида, находяще- 
гося в части системы Ф. При определении объе- 
мной активности /-ого радионуклида на поверхно- 
сти породы йф+з, контактирующей с частью систе- 
мы Ф= 1,2,3 полагается, что доля поверхности, с 
которой контактирует часть системы Ф, пропор- 
циональна ее насыщенности: 


й 


Ф 


А‘ 

(і-»ога ф _ 3 ’ 


(ф = 4-6). 


Количество воды в газовой фазе описывается с 
помощью объемной концентрации С/: 
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с ,_м\_ М\ 

3 Ѵ 3 тѴ8 3 ' 


где Мф - масса /-ого компонента, находящегося в 
части системы Ф. 

Вертикальная фильтрация жидкой фазы опре- 
деляется капиллярными и гравитационными сила- 
ми. Значение скорости фильтрации (/, жидкой фа- 
зы в ненасыщенной зоне рассчитывается с помо- 
щью уравнения Клюта-Ричардса [3]: 


где к - абсолютная проницаемость;^,) - функция 
фазовой проницаемости жидкой фазы; /л - вяз- 
кость воды; § - ускорение свободного падения; 
*Р - потенциал капиллярного давления. В настоя- 
щей работе для определения зависимости фазовой 
проницаемости /(5,) от насыщенности 5, использу- 
ется степенная функция. Кроме этого полагается, 
что движение сплошной водной фазы может про- 
исходить при условии, что ее насыщенность пре- 
вышает некоторое критическое значение 5,’ кото- 
рое называется остаточной насыщенностью [4] : 


Ж) = 


Ч-*; ѵ 

1-5,* 

О, 


5,>5; 

5,<5,* 


где у- параметр модели. В работе используется ли- 
нейная зависимость потенциала капиллярного да- 
вления от свободной влажности [3, 5]: 


Ч І = Ь-і 


ч-^; л 

т - 5,* 


где Ъ, Щ. - параметры потенциала капиллярного да- 
вления. 


Скорость вертикальной фильтрации газовой 
фазы определяется из закона сохранения объе- 
ма, записанного с учетом несжимаемости жидкой и 
твердой фаз: 


ад 


-С 1+ { 


4 

р‘ 


Т ' 3 КТ ^ т8 ъ ѲТ 8Р Л С 

М'Р + Т ді Р ді ) 


где Т - температура; Р— давление; р‘ - плотность 
/-ого компонента; - плотность потока /-ого 
компонента из части Ф в часть Ф; ! — время; 
К - универсальная газовая постоянная; М' - мо- 
лярная масса /-ого компонента. 

Плотность потока воды , между жидкой и га- 
зовой фазами в результате испарения (конденса- 
ции) пропорциональна отклонению парциального 
давления водяного пара в газовой фазе Р 1 от равно- 
весного значения Р'„( Т): 

Л п =у ] п 8 ]Ъ (Р^Т)-Р'), 


где (5 ф д, - удельная площадь контакта частей систе- 
мы Ф и Ф. Используя уравнения состояния Менде- 
леева-Клайперона для водяного пара, можно найти 
величину парциального давления Р'\ 


ы _ С\КТ 

м 1 

Плотность потока воды /) 2 между жидкой фазой и 
льдом рассчитывается из баланса тепловой энергии. 

Молекулярная диффузия радионуклидов в жид- 
кой фазе и водяного пара в газовой фазе описыва- 
ется на основе закона Фика. Диффузионный поток 
остается пропорциональным градиенту концентра- 
ции. Для учета замедления диффузии вследствие 
изогнутости пор и наличия в системе других фаз в 
законе Фика используется эффективный коэффи- 
циент диффузии /-ого компонента в части системы 
Ф Офм, связанный с коэффициентом молекуляр- 
ной диффузии в фазе Ф соотношением [6]: 

_ Д>5 ф 

и ф ,Эф ’ 


где х ~ извилистость. 

Плотность потока /| 4 /-ого радионуклида между 
жидкой фазой и поверхностью вмещающей поро- 
ды в результате адсорбции/десорбции рассчитыва- 
ется в рамках приближения линейной кинетики, 
согласно которой /; 4 пропорциональна отклоне- 
нию от равновесного значения, определяемого с 
помощью изотермы Генри: 

Л,4 = 71,4^1,4 ( а 1 “ ^1,4^4 )• 


где К\ А - коэффициент распределения для /-ого 
компонента между жидкой фазой и поверхностью 
вмещающей породы; у\ 4 - постоянная скорости 
массообмена между жидкой фазой и поверхностью 
вмещающей породы для /-ого компонента. 

Изменение насыщенностей частей системы с 1 
по 3, будет сопровождаться переходом части ГРФ 
от одной фазы к другой и изменением объемных 
активностей радионуклидов в них. Для описания 
этого процесса вводятся эффективные массопото- 
ки Л А5 и /) б . Значение потока /) 5 пропорционально 
изменению насыщенности льда: 
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а /Д - изменению насыщенности газовой фазы: 
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где Ѳ(х) - ступенчатая функция (Ѳ(х)=0 при х>0 и 
Ѳ(х)=0 прих<0). 

Расчет динамики температурного поля основы- 
вается на предположении о локальном тепловом 
равновесии между всеми частями системы. Объе- 
мная теплоемкость среды определяется ее составом: 

3 

с = с 0 (1-т) + т^8 ф с ф , 

Ф=1 

где с ф - объемная теплоемкость части системы Ф. 
Теплопроводность Д системы определяется уравне- 
нием Максвелла [7]: 
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Р 


ЗА) -2 т Д.+2ЗД 


3 + т 


РоГЕ$ Ф Р Ф - 1 


где - теплопроводность части системы Ф. 

Энерговыделение в системе происходит в резуль- 
тате радиоактивного распада и фазовых переходов. 
Мгновенное значение удельной мощности \Ѵ' энер- 
говыделения в результате распада /-ого радиоактив- 
ного радионуклида рассчитывается по формуле: 

IV' =^Ѵа;„, 

Ф 


где \ѵ‘ - энергетическая постоянная для /-ого ком- 
понента. Удельная мощность энерговыделения \Ѵ 1 
в результате испарения/конденсации и кристалли- 
зации/плавления воды определяется формулой: 

\У'=Ь і'-і У, 1 ,, 

пл 1,2 йен 1,3 ? 


где Ь ш и Ь исп - скрытые теплоты кристаллиза- 
ции/плавления и испарения/конденсации воды. 

Таким образом, изменение насыщенностей ча- 
стей системы определяется следующей системой 
дифференциальных уравнений: 
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Изменение объемной активности /-ого радиону- 
клида рассчитывается из закона сохранения массы: 

г л л 
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где 2' - постоянная распада /-ого радионуклида. 
Динамика концентрации воды в газовой фазе 
определяется уравнением: 
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Уравнение для расчета динамики температур- 
ного поля системы, записанное с учетом конвек- 
тивного теплообмена, теплопроводности и энерго- 
выделения в результате радиоактивного распада и 
фазовых переходов, имеет вид: 
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Для прогнозирования миграции радионуклидов 
в поверхностном слое почвы необходимо учитывать 
внешнее воздействие со стороны окружающей сре- 
ды. К основным факторам внешнего воздействия 


относят: изменение атмосферного давления, темпе- 
ратуры и влажности воздуха, атмосферные осадки. 
Естественные факторы (давление, температура, 
влажность) подвержены периодическим колеба- 
ниям (суточные, годовые циклы), а также случай- 
ным изменениям в соответствии с меняющимися 
погодными условиями. В настоящей модели для 
определения атмосферного давления, температуры 
поверхности земли и влажности атмосферного воз- 
духа используются периодические граничные усло- 
вия, описывающие суточные и сезонные колебания 
этих величин и позволяющие учесть внешнее воз- 
действие на рассматриваемую систему. Суммарное 
количество осадков, выпавшее на поверхность поч- 
вы за один месяц, соответствовало данным, пред- 
ставленным в работе [8]. При этом интенсивность, 
продолжительность и время начала осадков задава- 
лись произвольным образом и соответствовало 
средней полосе России. В модели учитывается уве- 
личение инфильтрации воды с поверхности весной 
в результате таяния снега. 

На основе численной реализации представ- 
ленной математической модели было разработано 
проблемно-ориентированное программное обеспе- 
чение, позволяющее моделировать вертикальную 
миграцию радионуклидов в поверхностном слое 
почвы. С помощью программного обеспечения про- 
ведено моделирование миграции Зг" и Ск 137 с поверх- 
ности почвы. Исследовались динамика температур- 
ного поля с учетом радиогенного энерговыделения, 
влияние различных физико-химических процессов 
(инфильтрации воды с поверхности, взаимодействие 
с вмещающей породой) на вертикальную миграцию 
радионуклидов. В начальный момент времени ра- 
дионуклиды находились в адсорбированном состоя- 
нии на поверхности почвы, в слое толщиной 5 мм. 
Суммарная активность радионуклидов на единицу 
площади составляла 0,5 мКи/м 2 . Коэффициенты ра- 
спределения между жидкой фазой и поверхностью 
породы - А) 5 4=1,35 , А) с )= 40,5. Пористость почвы 
м=0,4. Остальные параметры моделирования опре- 
делялись по данным [8, 9]. 

На динамику температурного поля основное 
влияние оказывают годовые и суточные изменения 
температуры поверхности почвы. Суточные и годо- 
вые колебания температуры поверхности приводят к 
формированию температурных волн, распространяю- 
щихся вглубь почвы. Прохождение температурной 
волны сопровождается колебаниями температуры 
грунта. По мере увеличения глубины температурные 
волны затухают. Интенсивность затухания тем боль- 
ше (соответственно глубина проникновения волн тем 
меньше), чем меньше период колебаний температуры 
на поверхности. Суточные колебания температуры 
приводят к изменению температуры только в припо- 
верхностном слое грунта. Амплитуда суточных коле- 
баний температуры уменьшается с глубиной и стано- 
вится пренебрежимо мала на глубине 30 см (рис. 1). 
Амплитуда годовых колебаний также уменьшается с 
глубиной. Поскольку период годовых колебаний 
больше периода суточных колебаний, то глубина про- 
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никновения температурных волн больше и составляет 
1,5 м (рис. 2). Полученные результаты моделирования 
динамики температурного поля в поверхностном слое 
почвы согласуются с экспериментальными наблюде- 
ниями [8] . Моделирование с различной начальной ак- 
тивностью радионуклидов показало, что радиогенное 
тепловыделение приводит к заметному повышению 
температуры в поверхностном слое почвы лишь при 
значительных активностях (более ІО 3 Ки/м 2 ). 



Рис. 1. Изменение температуры на различной глубине в те- 
чение двух суток ( июль): 1) на поверхности, 2) на глу- 
бине 5, 3 ) 10, 4) 15, 5) 20 см 



Время, месяц 

Рис. 2. Г одовые колебания температуры на различной глу- 
бине (в 12 ч дня): 1) на поверхности, 2) 20, 3) 50, 4) 
100 см 



Рис. 3. Миграция стронция с конвективным потоком жидкой 
фазы: 1 - 4 ~ насыщенность водной фазы, 5-8 ~ объе- 
мная активность стронция в водной фазе; 1, 5) до вы- 
падения осадков; 2, 6) через 3 ч, 3, 7) через 40 ч, 4, 
8) через 70 ч после начала выпадения осадков 



Рис. 4. Объемная активность стронция в жидкой фазе: 
1) 1 год, 2) 3 года, 3) 5 лет, 4) 7 лет, 5) 10 лет 



Рис. 5. Объемная активность цезия в жидкой фазе: 1) 1 год, 
2) 10 лет, 3) 20 лет, 4) 30 лет 

Вертикальная миграция радионуклидов происхо- 
дит в результате молекулярной диффузии и конвек- 
тивного движения жидкой фазы. Отсутствие жидкой 
фазы в пористой среде в результате ее замерзания 
или испарения приводит к прекращению миграции 
радионуклидов. Вследствие малой скорости диффу- 
зии радионуклидов в поверхностном слое почвы ос- 
новным механизмом вертикальной миграции явля- 
ется массоперенос с нисходящим потоком жидкой 
фазы. На рис. 3 приведены результаты моделирова- 
ния миграции стронция в результате инфильтрации 
воды с поверхности почвы после выпадения осадков. 
До выпадения осадков насыщенность воды в почве 
определяется капиллярными силами и составляет 
10 % от объема пор (кривая 1). Объемная активность 
стронция в жидкой фазе имеет форму вала, сформи- 
ровавшегося в результате миграции с водой в течение 
трех лет после попадания стронция на поверхность 
почвы (кривая 5). После выпадения осадков, в почве 
формируется область с повышенной насыщенно- 
стью водной фазы. Наибольшая насыщенность 
определяется интенсивностью осадков, а размеры 
области с повышенной насыщенностью - суммар- 
ным объемом выпавших осадков (кривая 2). Увели- 
чение насыщенности водной фазы в области залега- 
ния радионуклидов вследствие десорбции приводит 
к увеличению их концентрации в жидкости (кри- 
вая 6). При движении водной фазы под действием 
гравитационных сил в результате сорбции происхо- 
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дит отставание фронта активности (кривые 7, 8) от 
фронта насыщенности (кривые 3, 4). Чем больше ко- 
эффициент распределения радионуклида между 
жидкой и твердой фазами, тем меньше его концен- 
трация в жидкой фазе и соответственно тем медлен- 
нее происходит его вертикальная миграция. На 
рис. 4 и 5 представлены вертикальные распределе- 
ния объемной активности стронция и цезия в раз- 
личные моменты времени. Большой коэффициент 
распределения стронция между жидкой фазой и по- 
родой приводит к замедлению его скорости мигра- 
ции. Таким образом, скорость вертикальной мигра- 
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